{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Racial data vs. Congressional districts\n",
    "\n",
    "We are now awash with data from different sources, but pulling it all together to gain insights can be difficult for many reasons.  In this notebook we show how to combine data of very different types to show previously hidden relationships:\n",
    "\n",
    "* **\"Big data\"**: 300 million points indicating the location and racial or ethnic category of each resident of the USA in the 2010 census.  See the [datashader census notebook](https://anaconda.org/jbednar/census) for a detailed analysis.  Most tools would need to massively downsample this data before it could be displayed.\n",
    "* **Map data**: Image tiles from ArcGIS showing natural geographic boundaries.  Requires alignment and overlaying to match the census data.\n",
    "* **Geographic shapes**: 2015 Congressional districts for the USA, downloaded from census.gov.  Requires reprojection to match the coordinate system of the image tiles.\n",
    "\n",
    "Few if any tools can alone handle all of these data sources, but here we'll show how freely available Python packages can easily be combined to explore even large, complex datasets interactively in a web browser.  The resulting plots make it simple to explore how the racial distribution of the USA population corresponds to the geographic features of each region and how both of these are reflected in the shape of US Congressional districts.  For instance, here's an example of using this notebook to zoom in to Houston, revealing a very precisely [gerrymandered](https://en.wikipedia.org/wiki/Gerrymandering_in_the_United_States) [Hispanic district](https://green.house.gov/about/our-district):\n",
    "\n",
    "![Houston district 29](../assets/images/houston_district29.png)\n",
    "\n",
    "Here the US population is rendered using racial category using the key shown, with more intense colors indicating a higher population density in that pixel, and the geographic background being dimly visible where population density is low.  Racially integrated neighborhoods show up as intermediate or locally mixed colors, but most neighborhoods are quite segregated, and in this case the congressional district boundary shown clearly follows the borders of this segregation.\n",
    "\n",
    "If you run this notebook and zoom in on any urban region of interest, you can click on an area with a concentration of one racial or ethnic group to see for yourself if that district follows geographic features, state boundaries, the racial distribution, or some combination thereof.\n",
    "\n",
    "Numerous Python packages are required for this type of analysis to work, all coordinated using [conda](http://conda.pydata.org):\n",
    "\n",
    "* [Numba](http://numba.pydata.org): Compiles low-level numerical code written in Python into very fast machine code\n",
    "* [Dask](http://dask.pydata.org): Distributes these numba-based workloads across multiple processing cores in your machine\n",
    "* [Datashader](http://datashader.readthedocs.io): Using Numba and Dask, aggregates big datasets into a fixed-sized array suitable for display in the browser\n",
    "* [GeoViews](http://geo.holoviews.org/) (using [Cartopy](http://scitools.org.uk/cartopy)): Project longitude, latitude shapes into Web Mercator and create visible objects\n",
    "* [HoloViews](http://holoviews.org/): Flexibly combine each of the data sources into a just-in-time displayable, interactive plot\n",
    "* [Bokeh](http://bokeh.pydata.org/): Generate JavaScript-based interactive plot from HoloViews declarative specification\n",
    "\n",
    "Each package is maintained independently and focuses on doing one job really well, but they all combine seamlessly and with very little code to solve complex problems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import holoviews as hv\n",
    "from holoviews import opts\n",
    "import geoviews as gv\n",
    "import datashader as ds\n",
    "import dask.dataframe as dd\n",
    "from cartopy import crs\n",
    "\n",
    "from holoviews.operation.datashader import datashade\n",
    "\n",
    "hv.extension('bokeh', width=95)\n",
    "\n",
    "opts.defaults(\n",
    "    opts.Points(apply_ranges=False, ),\n",
    "    opts.RGB(width=1200, height=682, xaxis=None, yaxis=None, show_grid=False),\n",
    "    opts.Shape(fill_alpha=0, line_width=1.5, apply_ranges=False, tools=['tap']),\n",
    "    opts.WMTS(alpha=0.5)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, we'll load data from different sources and show it all overlaid together.  First, let's define a color key for racial/ethnic categories:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "color_key = {'w':'blue',  'b':'green', 'a':'red',   'h':'orange',   'o':'saddlebrown'}\n",
    "races     = {'w':'White', 'b':'Black', 'a':'Asian', 'h':'Hispanic', 'o':'Other'}\n",
    "\n",
    "color_points = hv.NdOverlay(\n",
    "    {races[k]: gv.Points([0,0], crs=crs.PlateCarree()).opts(color=v) for k, v in color_key.items()})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we'll load the 2010 US Census, with the location and race or ethnicity of every US resident as of that year (300 million data points), and define a plot using datashader to show this data with the given color key:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = dd.io.parquet.read_parquet('../data/census.snappy.parq')\n",
    "df = df.persist()\n",
    "census_points = hv.Points(df, kdims=['easting', 'northing'], vdims=['race'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can datashade and render these points, coloring the points by race:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_range, y_range = ((-13884029.0, -7453303.5), (2818291.5, 6335972.0)) # Continental USA\n",
    "shade_defaults = dict(x_range=x_range, y_range=y_range, x_sampling=10, y_sampling=10, width=1200, height=682,\n",
    "                      color_key=color_key, aggregator=ds.count_cat('race'),)\n",
    "shaded = datashade(census_points, **shade_defaults)\n",
    "shaded"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we'll load congressional districts from a publicly available [shapefile](https://www.census.gov/geo/maps-data/data/cbf/cbf_cds.html) and project them into Web Mercator format using GeoViews (which in turn calls Cartopy):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shape_path = '../data/cb_2015_us_cd114_5m.shp'\n",
    "districts = gv.Shape.from_shapefile(shape_path, crs=crs.PlateCarree())\n",
    "districts = gv.project(districts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we'll define some image tiles to use as a background, using any publicly available Web Mercator tile set:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tiles = gv.WMTS('https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{Z}/{Y}/{X}.jpg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each of these data sources can be visualized on their own (just type their name in a separate cell), but they can also easily be combined into a single overlaid plot to see the relationships:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opts.defaults(\n",
    "    opts.Polygons(fill_alpha=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shaded = datashade(census_points, **shade_defaults)\n",
    "tiles * shaded * color_points * districts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should now be able to interactively explore these three linked datasets, to see how they all relate to each other.  In a live notebook, this plot will support a variety of interactive features:\n",
    "\n",
    "* Pan/zoom: Select the \"wheel zoom\" tool at the left, and you can zoom in on any region of interest using your scroll wheel.  The shapes should update immediately, while the map tiles will update as soon as they are loaded from the external server, and the racial data will be updated once it has been rendered for the current viewport by datashader.  This behavior is the default for any HoloViews plot using a Bokeh backend.\n",
    "* Tapping: click on any region of the USA and the Congressional district for that region will be highlighted (and the rest dimmed).  This behavior was enabled for the shape outlines by specifying the \"tap\" tool in the options above.\n",
    "\n",
    "Most of these interactive features are also available in the static HTML copy visible at [anaconda.org](https://anaconda.org/jbednar/census-hv-dask), with the restriction that because there is no Python process running, the racial/population data will be limited to the resolution at which it was initially rendered, rather than being dynamically re-rendered to fit the current zoom level.  Thus in a static copy, the data will look pixelated, whereas in the live server you can zoom all the way down to individual datapoints (people) in each region."
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
